<<< Only Problem 1 and 2 will be graded >>>¶

Problem 1 (sound)¶

Denoising time with FFT (DFT)

In [3]:
!pip install praat-parselmouth
Collecting praat-parselmouth
  Downloading praat_parselmouth-0.4.5-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl.metadata (2.9 kB)
Requirement already satisfied: numpy>=1.7.0 in /usr/local/lib/python3.10/dist-packages (from praat-parselmouth) (1.26.4)
Downloading praat_parselmouth-0.4.5-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (10.7 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 10.7/10.7 MB 55.1 MB/s eta 0:00:00
Installing collected packages: praat-parselmouth
Successfully installed praat-parselmouth-0.4.5
In [4]:
import numpy as np
import pandas as pd
from scipy import signal,fftpack
import cv2
from skimage.io import imread
import matplotlib.pyplot as plt
import IPython.display as ipd
import os

import librosa
import parselmouth
In [6]:
sampling_rate = 32000
N=10001
Nf = 3 # Nf--> num freq
t= np.arange(N,dtype=float)
# pick rand period betwwen 10-2010 and convert to freq

# random period
Ts = np.random.rand(Nf)*2000+10
fs=1/Ts

# fs in sampling rate = 32000
fs_real = fs*sampling_rate

# pick rand Amp and phase
amp = np.random.rand(Nf)*200+ 100
phi = np.random.rand(Nf)*2*np.pi

# create clean signal
h = np.zeros(N)
for i in range(len(fs)):
    h += amp[i]*np.sin(2*np.pi*fs[i]*t + phi[i])

# signal with noise
h_w_noise = h + np.random.randn(N)*3*h + np.random.randn(N)*700
In [7]:
# TODO 1.1 : plot (1) clean signal and (2) noisy signal with label

plt.figure(figsize=(15, 5))
plt.plot(t, h, label='Clean Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.title('Clean Signal')
plt.show()

plt.figure(figsize=(15, 5))
plt.plot(t, h_w_noise, label='Noisy Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.title('Noisy Signal')
plt.show()
No description has been provided for this image
No description has been provided for this image
In [8]:
# TODO 1.2: plot magnitude of FFT of the noisy signal (freq sort form min--> max)
from scipy.fft import fft, fftfreq, ifft
freqs = fftfreq(N, 1 / sampling_rate)
fft_values = fft(h_w_noise)


sorted_indices = np.argsort(freqs)
sorted_freqs = freqs[sorted_indices]
sorted_magnitude = np.abs(fft_values[sorted_indices])

plt.figure(figsize=(15, 5))
plt.plot(sorted_freqs, sorted_magnitude)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('Magnitude of FFT of Noisy Signal')
plt.show()
No description has been provided for this image
In [17]:
# TODO 1.3 : cleaning the noisy signal using magnitude of FFT

threshold = np.max(np.abs(fft_values)) * 0.5  # Threshold (50% of max value)
fft_cleaned = np.where(np.abs(fft_values) > threshold, fft_values, 0)  # Keep significant frequencies
h_filtered = np.real(ifft(fft_cleaned))


plt.figure(figsize=(15, 5))
plt.plot(t, h_filtered, label='Filtered Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.title('Filtered Signal')
plt.show()
No description has been provided for this image
In [14]:
# TODO 1.4 : plot clean signal, noise signal and filtered signal (from your result in TODO 3.3) with label

plt.figure(figsize=(15, 5))
plt.subplot(3, 1, 1)
plt.plot(t, h, label='Clean Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.title('Clean Signal')

plt.subplot(3, 1, 2)
plt.plot(t, h_w_noise, label='Noisy Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.title('Noisy Signal')

plt.subplot(3, 1, 3)
plt.plot(t, h_filtered, label='Filtered Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.title('Filtered Signal')

plt.tight_layout()
No description has been provided for this image
In [15]:
# TODO 1.5 : export with IPython.display, listen to (1) original signal (2) signal with noise (3) signal after filtered

ipd.display(ipd.Audio(h, rate=sampling_rate))
ipd.display(ipd.Audio(h_w_noise, rate=sampling_rate))
ipd.display(ipd.Audio(h_filtered, rate=sampling_rate))
Your browser does not support the audio element.
Your browser does not support the audio element.
Your browser does not support the audio element.
In [ ]:
# TODO 1.6 : Write to explain and analyze the results

Problem 2 (image FFT)¶

Download a 1000 x 1000 image ("hamtaro.png") below

hamtaro

In [20]:
screen_shot = cv2.imread('/content/hamtaro.png',0)

plt.figure(figsize=(10, 10))
plt.imshow(screen_shot, cmap='gray')
plt.show()
No description has been provided for this image
In [21]:
# Apply FFT to the given image
F1 = fftpack.fft2((screen_shot).astype(float))
F2 = fftpack.fftshift(F1) # FFT center zeros freq
plt.figure(figsize=(10,10))
plt.imshow( (20*np.log10( 0.1 + np.abs(F2))).astype(int), cmap=plt.cm.gray)
plt.show()
No description has been provided for this image
In [22]:
# TODO 2.1 : Implement an ideal high-pass filter with a box size of 100x100 pixels on the given image
n = 100
# Get the number of rows and columns in the image
rows, cols = screen_shot.shape
crow, ccol = rows // 2 , cols // 2

# Create a mask with a box of 100x100 pixels set to 0 (high-pass filter)
mask = np.ones((rows, cols), np.uint8)
mask[crow - n // 2:crow + n // 2, ccol - n // 2:ccol + n // 2] = 0

# Apply the mask to the shifted FFT
F2_high_pass = F2 * mask

# Inverse FFT to get the filtered image
F1_high_pass = fftpack.ifftshift(F2_high_pass)
img_high_pass = np.abs(fftpack.ifft2(F1_high_pass))

plt.figure(figsize=(10, 10))
plt.imshow(img_high_pass, cmap='gray')
plt.title('High-Pass Filtered Image')
plt.show()
No description has been provided for this image
In [24]:
# TODO 2.2 : Implement an ideal low-pass filter with a box size of 100x100 pixels on the given image
# Create a mask with a box of 100x100 pixels set to 1 (low-pass filter)
mask = np.zeros((rows, cols), np.uint8)
mask[crow - n // 2:crow + n // 2, ccol - n // 2:ccol + n // 2] = 1

# Apply the mask to the shifted FFT
F2_low_pass = F2 * mask

# Inverse FFT to get the filtered image
F1_low_pass = fftpack.ifftshift(F2_low_pass)
img_low_pass = np.abs(fftpack.ifft2(F1_low_pass))

plt.figure(figsize=(10, 10))
plt.imshow(img_low_pass, cmap='gray')
plt.title('Low-Pass Filtered Image')
plt.show()
No description has been provided for this image

Problem 3¶

A digital signal can be generated from sampling of an analog signal using a periodic impulse-train. Explain how you can reconstruct an analog signal from a digital signal and aliasing problem does not occur when $f_s \leq 2f_{max} $ using frequency analysis.

where $f_s$ is the sampling frequency and $f_{max} $ is the maximum frequency of the analog signal

HINT : $ \mathscr{F} \left\{ \sum_{n=-\infty}^{\infty} \delta (t-n T_s) \right\} = \sum_{n=-\infty}^{\infty} \delta(\omega - n\omega_s)$ if $\omega_s = \frac{2\pi}{T_s} = 2\pi f_s$

Problem 4 : Aliasing¶

Problem 4.1¶

The following code generates two sine waves (x01_ts01 and x02_ts01) which are sampled in a range of t = 0,0.05 with sampling rate = 5000 Hz (f_samp_01). Study and write a report to analyze the results.

In [ ]:
t_st = 0
t_end = 0.05
f_01 = 200
f_02 = 2300

f_samp_01 = 5000

ts01 = np.linspace(t_st, t_end , int((t_end-t_st)*f_samp_01), endpoint=False)
x01_ts01 = np.sin(2*np.pi*f_01*ts01)
x02_ts01 = np.sin(2*np.pi*f_02*ts01)

plt.figure(figsize=(20, 5))
plt.plot(ts01, x01_ts01, 'go-', ts01, x02_ts01, 'r.-')
plt.show()
No description has been provided for this image

The sampling rate is reduced to 2500 Hz (f_samp_02). Study and write a report to compare the results.

In [ ]:
f_samp_02 = 2500
ts02 = np.linspace(t_st, t_end , int((t_end-t_st)*f_samp_02), endpoint=False)
x01_ts02 = np.sin(2*np.pi*f_01*ts02)
x02_ts02 = np.sin(2*np.pi*f_02*ts02)

plt.figure(figsize=(20, 5))
plt.plot(ts02, x01_ts02, 'go-', ts02, x02_ts02, 'r.-')
plt.show()
No description has been provided for this image

Ans.

Problem 4.2¶

The following code generate audio signals at different frequencies. Play the sound and write a report the analyse the results.

In [ ]:
t_st = 0
t_end = 5
f_01 = 50
f_02 = 22050 - f_01
f_03 = 22050 + f_01
f_samp_02 = 22050

ts02 = np.linspace(t_st, t_end , int((t_end-t_st)*f_samp_02), endpoint=False)

# CREATE SIGNAL WITH DIFFERENT FREQ

x01_ts02 = np.sin(2*np.pi*f_01*ts02)
x02_ts02 = np.sin(2*np.pi*f_02*ts02)
x03_ts02 = np.sin(2*np.pi*f_03*ts02)
In [ ]:
x02_ts02
In [ ]:
ipd.Audio(x01_ts02, rate=f_samp_02)
Out[ ]:
Your browser does not support the audio element.
In [ ]:
ipd.Audio(x02_ts02, rate=f_samp_02)
Out[ ]:
Your browser does not support the audio element.
In [ ]:
ipd.Audio(x03_ts02, rate=f_samp_02)
Out[ ]:
Your browser does not support the audio element.

No description has been provided for this image

from Imgflip Meme Generator

TODO : write report¶

Ans:

Problem 4.3¶

why many of audio file use sampling rate 44.1 kHz¶

Ans:

Problem 5¶

Download the 3 audio files and analyze all 3 signals with preliminary analysis. (HINT : Use a log scale for both frequency and magnitude.)"

  1. bass-guitar-single-note --> mixkit-bass-guitar-single-note-2331.wav

explain pattern of signal

In [ ]:
!wget https://raw.githubusercontent.com/Pataweepr/ComEngMath2_2023_resource/master/mixkit-adult-sneeze-2212.wav
!wget https://raw.githubusercontent.com/Pataweepr/ComEngMath2_2023_resource/master/mixkit-child-deep-breath-2237.wav
!wget https://raw.githubusercontent.com/Pataweepr/ComEngMath2_2023_resource/master/mixkit-bass-guitar-single-note-2331.wav
In [ ]: